Ring Exchange and Phase Separation in the two-dimensional Boson Hubbard model 
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We present Quantum Monte Carlo simulations of the soft-core bosonic Hubbard model with a ring 
exchange term K. For values of K which exceed roughly half the on-site repulsion U, the density 
is a non-monotonic function of the chemical potential, indicating that the system has a tendency 
to phase separate. This behavior is confirmed by an examination of the density-density structure 
factor at small momenta and real space images of the boson configurations. Adding a near-neighbor 
repulsion can compete with phase separation, but still does not give rise to a stable normal Bose 
liquid. 
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Introduction 

The most commonly studied Hamiltonians of quantum 
spins and bosons include kinetic energy in the form of two 
site exchange. However, it was realized early on that, for 
lighter atoms, larger numbers of particles can permute 
their positions, leading to ring exchange termsi Inter- 
est in such models subsequently developed for a num- 
ber of reasons, most recently as a result of the possi- 
bility of the existence of a normal bose liquid. Such a 
phase which might reconcile apparently contradictory ob- 
servations in cuprate superconductors in the pseudogap 
regime where there is evidence of local pairing without 
superconductivity^*^ Indeed, the magnitude of ring ex- 
change terms in high-T c materials could be a significant 
fraction of the two site exchange 

Over the past several years, considerable numerical ef- 
fort has been expended in looking at the low temperature 
phases of Hamiltonians including ring exchange .LSaLiil 
Interestingly, while the normal bose liquid has remained 
somewhat elusive, a rich panoply of other behavior has 
been observed, including striped order, Noel antiferro- 
magnetism, valence bond solids, and phase separation. 
The nature of the phase diagram and of the phase tran- 
sitions between the different ordered phases remains an 
area of active inquiry^! Very recently it was proposed 12 
that ring exchange terms might be engineered into cold 
atomic gases, opening up the possibility that the effects 
of such terms could be studied in a situation where they 
can be tuned systematically, and hence the theoretical 
predictions checked experimentally. 

In this paper we will extend our previous work& on 
phase separation in the boson Hubbard model with ring 
exchange. We focus on two dimensions, since that is the 
lowest dimension allowing ring exchange processes, and 
it is the effective dimension of the CuC>2 sheets of the 
cuprate superconductors which have motivated much of 
the recent study of Bose liquid phases. Our key con- 
clusion is that when the ring exchange energy scale ex- 
ceeds approximately one half the on-site repulsion, the 



ground state is thcrmodynamically unstable to phase sep- 
aration. Adding a near-neighbor repulsion can prevent 
phase separation, but the ground state always has ei- 
ther nonzero checkerboard or superfluid order. Because 
we use a break-up of the Hamiltonian in constructing 
the path integral which has not been much employed be- 
fore, we also describe some of the technical details of the 
simulation, including the effectiveness of different local 
Quantum Monte Carlo (QMC) moves. 

The most simple boson Hubbard modelM is: 



H = -t (a]oj + a\a rj ) + U ' 



\ni{fhi-l) (1) 



The operators a}- , a,j create (destroy) a boson on site 



j, and obey commutation rules [o^,at] = Sij. The 

number operator is hj = OjCij. The hopping parameter 
t measures the kinetic energy and U the strength of the 
on-site repulsion. The sum (ij) is over near neighbors on 
a square lattice. In the limit U — > oo this model maps 
onto the spin-1/2 XY model, with the magnetization 
playing the role of the bosonic density. By now, the 
T — phase diagram of the boson Hubbard Hamiltonian 
is well known ,iiiiiiSiiSiilii2ii£ Away from commensurate 
filling, the ground state is superfluid. At commensurate 
fillings, and for sufficiently large U, a Mott insulator 
forms in which each site is occupied by an integer 
number of particles. 

Our goal is to consider the effect of a ring exchange 
term: 

K, — — K Y^(aia\a\a,i + a\a2a^a\) (2) 
W 

The ring exchange term acts on four site plaquettes {p}, 
destroying two particles which lie along one diagonal, and 
creating them on the other. Because of the minus sign 
in HJ), the energy decreases when the particles are ex- 
changed. The basic qualitative picture behind the sug- 
gestion that the ring exchange term might give rise to a 
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FIG. 1: The ring exchange term acts on each square plaquette 
of our 2d lattice (a), and allows the exchange of two particles 
from one diagonal (b) to the other (c). 



normal Bose liquid is that if one starts with a superfluid 
phase, K introduces local vortices in which two particles 
jump in opposite directions. These local twists, if suf- 
ficiently favored energetically, could compete with, and 
ultimately prevent, a coherent long range superflow of 
particles. At the same time, such a hopping process in- 
creases (local) quantum fluctuations, and would not be 
expected to make the system incompressible (a Mott in- 
sulator). Therefore one might have a normal Bose liquid. 

As we shall see, we find that this ring exchange term 
causes phase separation rather than a normal Bose liq- 
uid. It is natural to attempt to resurrect the liquid by 
adding a longer range repulsion which acts against parti- 
cle clustering. The most simple form is a near-neighbor 
repulsion: 



V = V 



2^ n ^ n 3 



(3) 



The effects of near neighbor repulsion in the absence of 
ring exchange have been well-studied in the context of the 
boson Hubbard Hamiltonian^S They are perhaps most 
simply ascertained by considering the hard-core limit and 
the equivalent spin model. As commented earlier, the bo- 
son kinetic energy corresponds to an XY spin coupling. 
V corresponds to a Z coupling. Thus when V < 2t the 
spin model is in the XY universality class (a boson su- 
perfluid) and for V > 2t the spin model is in the Ising 
universality class (a checkerboard charge density boson 
insulator) . A considerable amount of work has also been 
done on looking for supersolid phases which combine su- 
perfluid and charge density wave (CDW) orderi 21 ; 22 



QMC: The Partition Function as a Path Integral 

Our Quantum Monte Carlo simulations are performed 
using a World Line algorithm with a decomposition in- 
volving four site matrix elements. We discuss the algo- 
rithm in detail here because this decomposition, which is 
similar but not identical to reference^, is used consider- 
ably less frequently than the more conventional approach 
which involves two site matrix elements^* In fact, to our 



knowledge, it has not been employed before for Hamilto- 
nians including ring exchange, although stochastic series 
expansion approaches 2 ^ are formally very similar. A spe- 
cific issue we will address is the effect of different local 
moves on the measurement of the energy. 

We begin by representing the partition function as a 
path integral, in which e~" is the imaginary time evo- 
lution operator, < r < (3. The goal is to take the 
exponential of the full Hamiltonian, which cannot be 
computed, and express it in terms of exponentials whose 
numerical values can be written down. In order to ac- 
complish this, we write the full imaginary time evolution 
operator as a product of infinitesimal evolution operators 
over short imaginary times At: 



Z = Tr 



-m - 



Tr (e- ArH ) M 



Here M is chosen to be a sufficiently large integer so as 
to make (At) 2 times any two of the energy scales in H 
small, as discussed further below. 

We then divide the sum over all plaquette operators 
into four classes depicted by the different shadings (colors 
on-line) in Fig. EI That is, the full Hamiltonian is written 
as 



where 



Tt = Tii + Tin + + H.4 



Ti n = ^2 ftp 



(5) 



(6) 



groups together all plaquettes of a given shading. In 
order to avoid overcounting the kinetic energy, on-site, 
and near neighbor interaction terms in the Hamiltonian 
we define, 



rip — K^p 



\(% + Vp) + \u p (7) 



with 



% = -*(ap 1 a p2 + a p3 a p4 + a pl a p3 + a p2 a p4 + h.c.) 
]Cp = -K (a p ia p2 a p3 a P 4 + h.c.) 

Up = U[h p i(n p i - 1) + h p2 (h p2 - 1) 

+h p3 (h p3 - 1) + n pi {n P 4 - I)] 
V p = V(n pl hp 2 + hpzhpi + n p ih p3 + h p2 h pi ). (8) 

It is important to emphasize that while plaquette oper- 
ators acting on neighboring plaquettes do not commute, 
all plaqette operators within the same shading do com- 
mute, since they do not "touch." This independence will 
enable us to compute their exponentials. 

Now we express each infinitesimal evolution operator 
as a product: 
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e -ATHi e -ATH2 e -ATH 3e -ATH4 



(9) 
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FIG. 2: The application of e~ propagates the boson 
occupations on the four spatial sites of the plaquette (a) in 
the imaginary time direction, and can be thought of geomet- 
rically as converting each spatial plaquette into a cube (b) 
whose bottom face contains the boson occupations before the 
action of the operator, and whose top face contains the boson 
occupations after the action of the operator. 



Substituting expression @ in expression (QJ, we get: 

The errors made in breaking up e~ Arn , Eq. [5J go as 
the commutator of the individual pieces, and therefore as 
(At) 2 . It might be thought that the accumulation of the 
(At) 2 errors across the M imaginary time slices might 
lead to an error linear in At. However, the expectation 
value of any Hermitian operator^S! has a Trotter error 
which is quadratic in At. 




(a) 



(b) 



FIG. 3: (a) Plaquettes with same shading (color online) 
are independent, (b) Introducing a dimension in imaginary 
time makes the plaquettes become cubes, each representing 
an infinitesimal evolution operator. Cubes with same shading 
don't touch each other, while cubes with different shading act 
at different times and therefore are decoupled. 

To complete the evaluation of the trace, we work in the 
basis of occupation numbers Inserting a complete 

set of states / = J2{ n } K n }) ({ n ll between each pair of 
exponentials, we get: 



k+1 



z =j2 n« 

{n k } k=l 

x ({n fc + 3 / 4 } 
x ({n fe+1 /2} 

fc+1/4 



e -AxW 4 || n fe+3/4^ 



x({n 



}|- 



-AtKi 



|{n fe+1 /2}) 
|| n fe+i/ 

IK}) 



-AtH 3 



e -ArH 2 || n fe+l/4^ 



(11) 



The superscripts k appearing in the sets {n fc } label the 
imaginary time, emphasizing that a different complete 
set of states has been introduced between each pair of 
exponentials. Since operators acting on plaquettes with 
the same shading commute, each matrix element in <|llf) 
can be written as a product of four site matrix elements 



IKI) 



=n<^ 

v{i) 



9+1/41 -ArH p 



\v q 



(12) 

where V q represents the state at imaginary time q (q is a 
multiple of j ) of the four sites belonging to plaquette p, 
and Y[ P (q) i s a product over all matrix elements between 
times q and q + 1/4. The entire evolution is represented 
by an interconnected set of the cubes of Fig.JSJa, as illus- 
trated in Fig. 12)3. 

Since each Hp conserves the sum of the boson numbers 
on the four spatial sites, the occupation numbers living 
on these cubes trace out continous 'world-lines' of parti- 
cles which wiggle and deform during their evolution. Fur- 
thermore, because we are evaluating the trace of e _/3w , 
the occupation numbers and associated world lines are 
periodic in imaginary time. We will discuss the sampling 
of the expression for Z, and the associated conservation 
laws, in more detail below. 

Summarizing, the partition function takes the simple 
form of a sum of products of four site matrix elements, 



Z = 



M 

Enn(p 



-AtH. 



\V k ) (13) 



The imaginary time index k runs from 1 to M with step 
i. Each matrix element in l|13l) can be computed by 
numerical diagonalization if we assume that occupation 
numbers are never greater than four. (This restriction is 
satisfied in all simulations here, since the average densi- 
ties studied are of order unityi2£) Each site has five pos- 
sible states and the number of different states for one 
plaquette is 5 4 = 625. To compute the exponential of Tip 
then requires the diagonalization of a 625 x 625 matrix. 



Sampling of the partition function 

In order to perform measurements, one needs to gener- 
ate sets of world lines with the Boltzmann weight given 
by the summand of Eq.^J This is done using a Metropo- 
lis algorithm. One starts with straight world lines, then 
suggests local deformations and accepts them with prob- 
ability p — min(l, 1Z), where 1Z is the ratio of the Boltz- 
mann weight after the move to that before. In the usual 
way, this satisfies detailed balance. To ensure ergodic- 
ity of the algorithm, different kinds of moves have to be 
included, like pull moves (Fig. 0Ji), "baby" pull moves 
(Fig.^Ja), twist moves (Fig.^t), and "baby" twist moves 
(Fig. 0Ji) . Only a few of the plaquette matrix elements 
change when one of these local moves is performed. The 



most complex is the twist move (Fig. |3J;) which involves 
10 matrix elements (only 5 cubes are shown in order to 
make the figure clearer). 




(a) (b) (c) (d) 



FIG. 4: To ensure ergodicity of the algorithm, four kinds of 
moves must be included: (a) pull moves, (b) "baby" pulls, (c) 
twist moves and (d) "baby" twist moves. 



2 



•3 

a. 



-4 - 





|3=0.25 All moves 


• • u=s 


0=0-25 No twist 


B£]u=s 


P=8 All moves 


■ ■ u=s 


P=8 No twist 


[ill 


2 


4 



K 

FIG. 6: Energy per particle for interacting bosons as a func- 
tion of the ring exchange strength K with and without in- 
cluding twist moves. 



These different kinds of moves are required in order 
to make the algorithm ergodic. As an illustration, we 
show in Fig. the energy per particle of free bosons as a 
function of the inverse temperature /3, with various types 
of moves suppressed. We can see that the correct energy 
E(/3) is obtained only when both pull and "baby" pull 
moves are included. 




FIG. 5: Energy per particle for free bosons as a function of 
the inverse temperature /3 with and without including pull 
and baby pull moves. 

It is notable in Fig. 5 that the correct energy is ob- 
tained even without the twist moves. It might be thought 
that this is because K = U = and that twist moves 
will become important when turning on interactions, es- 
pecially with ring exchange processes. However this is 
not the case, as can be seen in Fig. |SJ We believe this to 
be an advantage of the four site decomposition over the 
two site one. The four site decomposition already ana- 
lytically includes twists in the evaluation of the matrix 
elements of the exponentials of the H n . 



While including only pull and "baby" pull moves seems 
to get the energy right on the lattice sizes studied above, 
our simulations are done with all moves included. These 
local moves all take the form of local distortions of ex- 
isting world-lines and hence cannot change p. Thus our 
code works in the canonical ensemble. It simulates what- 
ever value of p with which we initialize the lattice. If 
desired, we can make contact with the grand-canonical 
ensemble, by computing the chemical potential as the en- 
ergy cost, in the ground state, of adding a boson to the 
system, p(N b ) — E(N b + 1) - E(N b ). Actually, canonical 
simulations have certain advantages in exploring phase 
separation, as we shall see, and as has previously been 
described^ 

Our simulations also are confined to the zero-winding 
sector of the space of all possible paths. This has im- 
plications for how we measure the supcrfluid density p s . 
Consider the well-known relation^ 



Ps 



(W 2 ) 
2dtp 



(14) 



where d is the dimension of the system, and (W 2 ) the 
mean square of the winding number operator which mea- 
sures the net flow of particles off the right side of the 
lattice and over to the left side. Our local moves cannot 
change the winding, as is evident by comparing the con- 
figurations of Fig. Hence, using Eq. (|14|l our algorithm 
would systematically give a zero value for the superfluid 
density, if we begin in a W = configuration. 

However, we can still measure p s . The procedure is 
as follows: We define the pseudo-current operator (using 
the notation of Fig. for the site labels) 

3x(k)= ^ {^cd + n c& - n c2 - n c4 ) 

ce(Kk+i) (15) 

3y{k) = 2^ (n C 7 + n cS - n c3 - n c4 ) 
ce(fc,fc-t-i) 
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FIG. 7: Example of world line with W = 1 (a) and W = 
(b) winding numbers in the case of a one dimensional lat- 
tice. Configurations with non-zero winding number cannot 
be reached with local moves. 



where the sum X) c e(fc fc+i) 1S over au cubes of Fig. [2b 
between times k and k + 1. These quantities measure 
the number of particles which jump in the positive x and 
y directions between imaginary times k and k + 1. The 
winding number operators in both x and y directions are 
then given by 



W x 



M 



N x 



k=l 
M 



Wy = 



^ v k=i 



(16) 



where N x and N y are respectively the number of sites 
in the x and y directions. We measure also the auto- 
correlation function of the pseudo-current: 



(17) 



J,(k) = J2\Uk + k%(k'] 

k' 

Jy(k)=J20y(k + k')jy(k / ) 



and compute its Fourier transform via 



(18) 



where u> n = -jjn. 



One can show that these pseudocurrent operators are 
related to the winding by, 



Jy(0) 

N 2 
y 



(19) 



As already stated, the measurement of the winding num- 
ber, and hence the Fourier transform at ui n = 0, will 
give a zero value, but we can avoid the problem by con- 
sidering the finite frequency response. In fact, one can 
showia that the correct, nonzero, winding number and p s 
are given by taking the zero frequency limit of J~ x {u> n ), 
as illustrated in Fig. |H1 
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FIG. 8: Restricting the winding number to a zero value intro- 
duces a discontinuity in the Fourier transform of the pseudo- 
current auto-correlation function. Performing an extrapola- 
tion gives the correct value of the superfluid density. 



Results: Hard-core limit, V — 

We can work in the hard-core limit by running a soft- 
core code at large U, or more directly by forbidding the 
acceptance of multiple occupancies in our Monte Carlo 
moves. This also simplifies the computation of the matrix 
elements in 1)130 since we can restrict to single occupancy 
states (and hence diagonalize only a 16 x 16 matrix). Such 
a hard-core model is equivalent to replacing the standard 
bosonic commutation rules by the hard-core rules, 



{a,, at} = 
{a,i,4} 





= 
= 1 



'% 5 ^J] 



J 



Vi^j 
-- f j 
-0 Vi^j 



(20) 



where {A, B) = AB + BA is an anti-commutator. The 
usual bosonic commutators between operators acting on 
different sites ensure that the wave function is symmetric 
in the exchange of particles, as required for bosons. 

Our simulations are done mainly on 16 x 16 lattices. 
We always take t = 1 for the hopping parameter, which 
thereby sets the energy scale. We will also set the near 
neighbor interaction V = until the final results section. 

In order to describe the different phases of the sys- 
tem, in addition to the superfluid density p s , we mea- 
sure the Fourier transforms of the density-density and 
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plaquette-plaquette order correlation functions, S(tt,tt), 
P(tt, 0) and P(0, tt) denned by, 



S\k x , kyj 



P[k x: ky} 



1 



(N x N y ) ^, 
1 

K 2 (N x N y y^, 



y^(rtf fi? + ,?,)e tk ' r 



(21) 



^(/C p /C p ,)e^-^-V)(22) 



Fig. El shows these quantities, and the superfiuid den- 
sity, as a function of if at half filling. Increasing K makes 
the superfiuid density decrease. p s vanishes completely 
around K = 8. At the same time, the plaquette struc- 
ture factors P(ir, 0) and P(0, tt) start to grow, indicating 
the formation of stripes by the ring exchange. This is the 
VBS (Valence Bond Solid) phase, previously discussed in 
quantum spin systems. Stripe formation is also observ- 
able by considering the correlations in the hopping from 
plaquette to plaquette, rather than the ring exchange. 
For higher values of K, the VBS phase becomes unstable 
and is replaced by a CDW as shown by the growth of 

Both the VBS and CDW phases are gapped. That is, 
the energy needed to add a particle to the lattice (the 
chemical potential) jumps abruptly across p — ^. We 
define G = p+ - /Lt_ with /i+ = E(N b = f + l) - £(f ) 
and ^_ = E(§) - B(f - l). G = in the superfiuid 
phase K < 8, and is non-zero for K > 8. The VBS 
and CDW phases are insulators. Fig. HOIshows G as one 
increases K out of the superfiuid phase at half-filling. 
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FIG. 9: The superfiuid density and the order parameters for 
the hard-core case at half filling for a 16 x 16 lattice. Increasing 
K destroys the superfiuid density. Then a solid order appears 
corresponding to a VBS phase. This phase becomes unstable 
if increasing K more, and yields to a CDW phase. 

We now consider what happens when the system is 
away from half-filling. Fig. ^] shows the energy as a 
function of density. There are several interesting fea- 
tures. First, the kink in the energy right at half-filling is 




FIG. 10: The gap in the chemical potential at half filling, 
as a function of K for a 16 x 16 Lattice. The gap starts to 
appear at if — 8, coincident with the vanishing of p s , and 
grows as K is increased, showing the presence of insulating 
phases. 




FIG. 11: The energy per site as a function of the density p, 
for if = 10 on a 8 x 8 Lattice. 



associated with a nonzero gap. Second, there is a region 
of densities for which the energy is concave down. As 
a consequence, the usual Maxwell construction indicates 
that it is energetically more favorable to have two sep- 
arate regions of different density than a single region of 
uniform density. The system phase separates. 

The energetic signature of phase separation is further 
illustrated in Fig.^J which shows the density of particles 
p as a function of the chemical potential p over a wide 
range of densities. The jump of the chemical potential at 
half filling is the gap associated with the solid VBS phase 
discussed above. The regions where the slope of the curve 
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FIG. 12: The density of particles as a function of the chemical 
potential, for K — 10 (8 x 8 Lattice). 



(proportional to the compressibility kt) is negative in- 
dicate that the system is thermodynamically unstable 22 
and undergoes phase separation. There are two such re- 
gions. The first is seen as soon as we go away from half 
filling, ie for the two fillings N/2±l immediately adjacent 
to half-filling, and indicates a first order transition from 
the VBS phase to a superfluid. (See also Fig. [Til ) The 
same negative compressibility near half-filling is obtained 
for values of K higher than 15, showing a first order tran- 
sition from the CDW phase to a superfluid. The second 
region of negative compressibility occurs below p c — 0.34. 

The spikes occuring at p sa 0.1 and p w 0.2 in Fig. 1121 
are not due to statistical fluctuations. Indeed, when 
phase separating, the system can form a stripe in its den- 
sity profile rather than an island, due to periodic bound- 
ary conditions. This produces a shift in the energy, and 
hence in the chemical potential. In our simulations, the 
system forms a stripe for p E [0.1,0.2], and an island 
when crossing the limits, generating the spikes. 

While analysis of the energy and chemical potential 
yield quantitative information about the regions of sta- 
bility, a more straightforward, qualitative, indication of 
clustering comes from simple real space images of the bo- 
son density during the course of a simulation. Whether 
clustering occurs or not, if the density is averaged over 
the whole simulation, we should observe a uniform den- 
sity distribution, since the probability of a state is in- 
dependant of any translation in the lattice. But if the 
density is averaged over only a few iterations, we observe 
a density profile reflecting the clustering fFig,IT3l. 

To corroborate further that this negative compressibil- 
ity is associated with phase separation, we define the ob- 
servable 51 to be the average of the density structure fac- 
tor over the smallest values of the wave vector: 



S{e x ,0) + S(e x ,e y ) + S(0,e y ) 



2tt 




FIG. 13: Typical QMC results for the average density distri- 
bution when the clustering occurs. For a density p = 12/256 
the system adopts a shape of an island, while for p — 48/256 
the system prefers to form a stripe, which, in the presence 
of periodic boudary conditions, is the geometry which maxi- 
mizes the number of occupied four site plaquettes for which 
the ring exchange energy is low. 



Q measures modulations of the density profile with 
wave lengths on the order of the size of the lattice. Fig. 1141 
shows the superfluid density and f2 as functions of the 
density of particles p, for a fixed value K = 10. The sym- 
metry of the curves reflects the particle-hole symmetry of 
the hard-core Hamiltonian. Starting from p = \ where 
Q and p s have a zero value, and which corresponds to 
the VBS phase, we see that doping the system increases 
the superfluid density rapidly. p s reaches its maximum 
for p ~ 0.4 (or p ~ 0.6 by symmetry). Then p s falls at 
the same time that grows and reaches its maximum 
for p ~ 0.2 (or p ~ 0.8). Fig. 1141 suggests that as soon 
as we go away from half-filling, there is a transition from 
the VBS phase first to a stable superfluid, but that sub- 
sequently phase separation sets in. The density p ~ 0.3 
where p s falls and increases rapidly matches quite well 
with the density at which k changes sign in Fig. 1121 
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FIG. 14: Superfluid density and SI as functions of the density 
of particles, for K = 10 on an 8 x 8 Lattice. 

The occurrence of phase separation indicated in the 
compressibility, f2, and in real space images, may be un- 
derstood with a simple physical picture. When K is suffi- 
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ciently large, the system increases the ring exchange pro- 
cesses in order to decrease its energy. But ring exchange 
is possible only if there are second neighbor particles. 
Thus in a dilute collection of particles, the ring exchange 
term has a tendency to act as an attractive potential. 
This is not true of the usual kinetic energy, since there a 
single boson can hop by itself without needing a partner. 

One can analyze this binding quantitatively for the 
case of two bosons by exact diagonalization (Lanczos). 
Fig. H5l shows the square root of the mean quadratic dis- 
tance (f 2 ) between two particles for different LxL lattice 
sizes in the ground state. The soft-core case is consid- 
ered with an on-site repulsion U = I2t. For low values of 
K the distance between the bosons grows linearly with 
L, as one would expect if the particles were distributed 
independently across the lattice. On the other hand, for 
high values of K, the distance between the particles is 
relatively independent of lattice size, suggesting the par- 
ticles stay together, and demonstrating the existence of 
a bound state. While the data shown are for U = 12t, 
the phenomenon is independent of U, and, in particular, 
is seen also in the hard-core case. The reason the on-site 
repulsion U does not successfully compete with the bind- 
ing is that the bound wavefunction favors the bosons on 
second-near-neighbor sites where U plays no role. This 
is reflected in the data of Fig. 1151 where we see that the 
distance between the paired bosons is still several lattice 
sites. 
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FIG. 15: The mean square distance y (f 2 ) between two 

particles as a function of K/t, in the soft-core case with U = 
12t at T = 0. 

Figs.ll2landll4lshowed our results for hard-core bosons 
at fixed K = 10 as the system is doped. By perform- 
ing further simulations for different values of K , we are 
able to generate the complete phase diagram of hard- 
core bosons (Fig. in the K — p plane. Gapped VBS 
and CDW phases exist only at half-filling and relatively 
large K. For weaker K and half-filling the system is 



superfluid. When the system is doped a small amount 
away from half- filling, only a superfluid phase is present. 
Upon further doping, the clustering region is reached. As 
K becomes large, clustering occurs closer and closer to 
half-filling. 

Melko et ali have studied the spin- 1/2 Hamiltonian, 

(ij) (ijkl) 

where 

Bjj = 2(SfSJ + SfS") ^ 
Viju = S i Sj S k S t + S i Sj S k 5, 

and S x , S v , S z , S + , S~ are the usual spin operators. 
This is the Heisenberg XY model with a four spin ex- 
change term and in a magnetic field. This model is ex- 
actly equivalent to our bosonic model in the hard-core 
limit. The magnetic field plays the role of the chemi- 
cal potential in the grand canonical ensemble. In the 
phase diagram obtained by Melko et ah, the superfluid, 
VBS and CDW phases are present and their boundaries 
match with our results. In addition the nature of the 
transitions also agree. Our simulations in the canonical 
ensemble exhibit clear phase separation indicating first 
order transitions from CDW insulator (Neel) to SF and 
also from SF to Mott (fully polarized, in the spin lan- 
guage). To see this, consider the phase diagram, Fig. 1 
in Ref. [7]. As discussed in that paper, the cuts along 
lines A and B show clear hysteresis indicating first or- 
der transitions. Our simulations in Fig. 1121 are done at 
constant K and would correspond most closely to a ver- 
tical cut, like A, in Ref. [7] but going from Neel to SF 
to fully polarized. This way two first order phase transi- 
tions would be crossed, as shown in Fig. 1121 here by the 
negative compressibility regions. Note also that the jump 
in the density when going from Neel to SF (e.g. cut B 
in [7]) is very small as shown in the magnetization his- 
togram, Fig. 3 of [7]. Recalling that p = m + 0.5, we see 
that the jump in p in the histogram is about 0.015, of 
the same order as we see in our Fig. 1121 

Thus, the hard-core bosonic model with ring exchange 
does not exhibit any compressible but non-superfluid, 
Bose-liquid, phase. In the next section we examine the 
possibility 2 that the relaxation of the hard-core con- 
straint could make the solid phases (VBS or CDW) evolve 
to a Bosc liquid. 

Results: Soft-core case, V — 

We now turn to the soft-core case, but still choose the 
intersite repulsion V — and half-filling. The superfluid 
density p s is shown as a function of K for U = St in 
Fig. El Ps first grows slightly when K increases. When 
K ~ it, p s starts to decrease. This decrease becomes 
rather rapid until, at K ~ 8i, p s levels off at about half 
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1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

FIG. 16: The phase diagram of hard-core bosons (16 x 16 
Lattice). The solid phases (VBS and CDW) exist only at 
half-filling. When doping the system, the phase first becomes 
superfluid, then the clustering occurs closer to p — ~ as K is 
increased. 

its K = value. This behavior is quite different from the 
hard-core case (inset) where p s decreases as soon as the 
ring exchange interaction is turned on, and vanishes at 
K ~ 8t. 
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FIG. 17: The superfluid density as a function of K for the 
soft-core case where U = 8t, at half-filling (16 x 16 Lattice). 
p s does not vanish when K is large, unlike the hard-core case. 

In order to understand the difference in the behavior 
of p s between the soft and hard core cases, we begin by 
examining O, the small momentum density-density struc- 
ture factor. Fig.ll8lshows VI as a function of K for U = At 
and U = 8t. In both cases, when K reaches a value of the 
order U/2, O grows sharply, showing, as in the hard-core 
case, the presence of a clustering tendency. Results for 
16 x 16 and 24x24 lattices are almost identical, indicating 
that the phase separation is not a finite size effect. 

The behavior of the superfluid density can then be ex- 
plained as follows: When K is weak, the ring exchange 
term supports the motion of the particles, and helps the 



kinetic term in delocalizing them. A flow can thus occur 
more easily and p s increases. When K is large, clustering 
occurs and long range flow across the lattice is somewhat 
inhibited. If, for example, the system forms a stripe, one 
might argue that roughly speaking the flow of particles 
is possible only in one direction, and p s should decrease 
by a factor of two, in agreement with Fig. 1171 




FIG. 18: n versus K for U = At and U = 8t and two different 
lattice sizes. The lattice is half-filled. We can see that Q 
increases sharply at K ~ U /2. 

As we have seen, at half-filling, in the hard-core case, 
superflow does not occur when K is large because of the 
formation of gapped VBS and checkerboard phases. For 
the soft-core case, no solid structure is established and 
the particles can circulate. Thus the soft-core model is 
always superfluid at half filling and undergoes a cluster- 
ing when K reaches a value of the order U/2. Fig. IT9l 
shows the phase diagram of soft-core bosons at half fill- 
ing. 




FIG. 19: The phase diagram of soft-core bosons (16 x 16 
Lattice). The phase is superfluid under the curve, while the 
system forms a cluster in the upper part. The curve is ap- 
proximately hyperbolic. 

Our phase diagram Fig. ^| is for half-filling, but the 
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main feature, a tendency to phase separation, is the same 
at other densities. As an extreme example, simulations 
done for a commensurate filling of one particle per site 
show that, like the superfluid phase, the Mott phase 
also collapses when the ring exchange processes are too 
strong, and is replaced by clustering (Fig. I20JI . 




FIG. 20: Density profiles for the soft-core case with U = 4 
and K = 8 for commensurate fillings, p = 1 (left) and p — 2 
(right). 



Results: Soft-core case, V =/= 

In the preceding two sections we have seen that a Bosc- 
liquid does not occur in either the hard- or soft-core 
Bose Hubbard models when ring exchange is included, 
in contrast to recent suggestions^ It is natural to con- 
sider whether longer range repulsion might prevent the 
system from collapsing. As a first step, we now include a 
repulsion between first nearest neighbors. Fig. 1211 shows 
the structure factor S(ir, ir) and £1 for U = 5t and differ- 
ent values of V, as a function of K at half-filling. Three 
different behaviors of the density correlations are observ- 
able: a CDW (density order at momentum (tt, ir)) when 
K is weak (well known in the model without ring ex- 
change processes) (21 at intermediate K a uniform phase 
where the density structure factor is small at all wave- 
lengths, and a regime of phase separation (il large) when 
K is big. As expected, V does suppress the tendency for 
clustering. K must be made larger for phase separation 
to occur when V is increased. The phases of Fig. 1211 arc 
also directly visible in real space snapshots of the densi- 
ties in the course of the simulations (Fig. l22*|) . 

Despite the absence of density order, the uniform phase 
(Fig. 122b ) is not a Bose liquid. Fig. [23 shows that the 
superfluid density is zero only for the CDW phase. As 
soon as the staggered density order is destroyed with in- 
creasing K, the superfluid density becomes nonzero. p s 
climbs rapidly and is roughly constant through the re- 
gion where both S(tt, it) and Q are small. It remains 
non-zero in the cluster phase, though it does drop by a 
factor of two when clumping begins. There is no Bose 
liquid phase in the model. (Of course, if the density is 
sufficiently small, then an island phase in which the bo- 
son cloud does not span the system will have vanishing 
p s . This situation also occurs for bosons in a confining 
potential, where the phase is nevertheless commonly re- 
ferred to as being superfluid. 
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FIG. 21: The CDW structure factor S{-k, tt) and Q for U = 5t 
and different values of the potential V between first nearest 
neighbors, at half filling (16 x 16 Lattice). V is seen to com- 
pete with phase separation, pushing the onset of the rise in Q 
out to larger and larger K. 
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FIG. 22: Density profiles at half filling for U = 5t, V = 5t, 
and different values of K: (a) K — t, CDW phase; (b) K = 5t, 
uniform density; and (c) K — 7t, phase separation. 
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FIG. 23: The superfluid density p 3 , the structure factor 
S(tv,tt), and f2 as functions of K for U = 5t and V — 5t 
(16 x 16 Lattice). The superfluid density vanishes only when 
the system is in the CDW phase. 
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Conclusions 

In this paper we have studied the effect of a ring ex- 
change term in the hard- and soft-core Bose Hubbard 
model, also including a near-neighbor repulsion, using 
Quantum Monte Carlo simulations in the canonical en- 
semble. It had been suggested that this term might 
lead to a normal Bose liquid phase, that is, one which 
is compressible but non supcrfluid. For the Hard-core 
case, we reproduced results obtained by Melko et aim in 
the grand canonical ensemble. However, working in the 
canonical ensemble enables us to capture and character- 
ize a hitherto unsuspected cluster phase of the model. 
No Bose liquid phase was observed. Thus the specula- 
tion of Paramekanti et aim, that the relaxation of the 
hard-core constraint might give rise to a Bose liquid does 
not appear to be borne out. The tendency to cluster- 
ing can be understood by analyzing the effect of K on a 
two boson system. We have shown numerically that K 
acts much like an attractive potential, and that when it 
is sufficiently large, the two boson form a bound state. 
This attraction leads to clustering in systems with larger 
boson numbers. 

We saw that, while a near neighbor repulsion competes 
with this clustering, it could not create a normal bose liq- 
uid. At weak K it promotes an alternate charge density 
wave insulator. At strong K, while V does help prevent 



clustering, p s remains nonzero. It might be interesting to 
study models with longer range interactions which will 
act together to prevent phase separation, but compete 
with each other by promoting phases at different order- 
ing momenta. Such models might, finally, realize the 
Bose liquid phase. 

To conclude, let us comment on the implications of 
our work for the phase diagram of the spin-1/2 quantum 
Heisenberg model with a ring exchange term. The kinetic 
energy term in the hard-core boson Hubbard model maps 
onto JJ2(SiS? + S?SJ) with exchange constant J = It. 
At the value U — At in our soft core model, double oc- 
cupancy is already very rare at half-filling, and hence we 
are almost in the hard-core limit. The value of the ring 
exchange energy scale required to drive phase separation 
for this U is K ks 2t, or in other words, K k, J . To re- 
produce the near-neighbor coupling of the z components 
of spin in the Heisenberg model we must include a near- 
neighbor repulsion in the bose-Hubbard model, a term 
which clearly would suppress phase separation. Thus we 
expect ring exchange to have the potential to drive phase 
separation in the Heisenberg model only for K consider- 
ably greater than J. 

We acknowledge support from the National Science 
Foundation under awards NSF DMR 0312261 and NSF 
INT 0124863, and useful input from R. Stones. 
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